Diagnosed Dissaggregation Try 2

Sequal to TryToDiagnoseDisagg.Rmd Which didn’t work the first time and kind of became a mess. And FluxInterceptSillyness.Rmd

24 March 2019

Start with core parts of the FluxInterceptSillyness solution

Then lets run dissaggregation on the flux profile.

Making a synthetic data set.

Data to make binned_EachSize_ald.csv come from Normalize_UVP_Flux.Rmd Which in turn pulls data generated by UVP2017.Rmd

options(readr.default_locale=readr::locale(tz="Mexico/General"))
source("UVP_2017_library.R")
#dataP2 <- bring_in_p2()
#twin01 <- make_twin_df_list(dataP2)
#ES01 <- twin01[[1]] %>% filter(depth <= 1000, profile == "stn_043")
#DS01 <- twin01[[2]] %>% filter(depth <= 1000, profile == "stn_043")
# # Modification, merge all stations and see what happens.
# ES01 <- twin01[[1]] #%>% filter(depth <= 1000)
# DS01 <- twin01[[2]] #%>% filter(depth <= 1000)

ES01 <- read_csv("dataOut/binned_EachSize.csv") %>% filter(depth <= 1000)
DS01 <- read_csv("dataOut/binned_DepthSummary.csv") %>% filter(depth <= 1000)
twinS <- list(ES01, DS01)

SimpleBins <- seq(from = 0, to = 1100, by = 100)

#debug(bin_depths)
binnedS <- bin_depths(twinS, bins = SimpleBins) %>% calc_psd

ESS <- binnedS[[1]]
DSS <- binnedS[[2]]

Different direction

Units of nparticles

# mySpec <- c(3.32405678509316, 1.85877856561554, 1.03359305917877, 0.564880885240349, 
# 0.309079523670959, 0.165424066706246, 0.0925270232480783, 0.0502491579183327, 
# 0.0275765808581977, 0.0147327962322043, 0.00847342228148116, 
# 0.00457079526918403, 0.00242657514481869, 0.00133840064588931, 
# 0.00073768231262027, 0.000406688319239427, 0.00021922892226164, 
# 0.000120831754575073, 6.6223357185171e-05, 3.58979192483917e-05, 
# 2.00799713729794e-05, 1.09115166969355e-05, 5.82645914172791e-06, 
# 3.28908275460212e-06, 1.57702715137e-06)

lb_vec <- c(0.102, 0.128, 0.161, 0.203, 0.256, 0.323, 0.406, 0.512, 0.645, 
0.813, 1.02, 1.29, 1.63, 2.05, 2.58, 3.25, 4.1, 5.16, 6.5, 8.19, 
10.3, 13, 16.4, 20.6, 26)

binsize_vec <- c(0.026, 0.033, 0.042, 0.053, 0.067, 0.083, 0.106, 0.133, 0.168, 
0.207, 0.27, 0.34, 0.42, 0.53, 0.67, 0.85, 1.06, 1.34, 1.69, 
2.11, 2.7, 3.4, 4.2, 5.4, 6)

# C_f = 4
#ag = 0.26 # as in rest of universe
C_f = 10.51 # for now
alpha = .52
gamma = .26
ag = alpha + gamma

C_f_global <- C_f
alpha_global <- alpha
gamma_global <- gamma
ag_global <- ag
# mySpec
# myDf <- tibble(lb_vec, mySpec, binsize_vec)

make_spectrum <- function(icp, psd, bins = lb_vec, binsizes = binsize_vec){
  C_n = exp(icp)
  nnp = (C_n * bins ^ psd) * (binsizes)
  nnp
}

testSpectrum <- make_spectrum( -2, -3.5)
testSpectrum

 [1] 1.038207e+01 5.952372e+00 3.394453e+00 1.903053e+00 1.068184e+00 5.865166e-01 3.364153e-01 1.874208e-01 1.055022e-01 5.781819e-02
[11] 3.409371e-02 1.887235e-02 1.028025e-02 5.814982e-03 3.287130e-03 1.858824e-03 1.027953e-03 5.810880e-04 3.266638e-04 1.816351e-04
[21] 1.041944e-04 5.808824e-05 3.182049e-05 1.841915e-05 9.060571e-06

Post for Stack overflow, at least in thory

library(tidyverse)
library(cowplot)

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************

particle_data <- ESS %>% select(-nparticles, -n_nparticles, -time, -binsize)

I have a dataset in which I have measured particles of different sizes at different depths. Each depth has several size bins associated with it lb: lower bound of particle size (mm) ub: upper bound of particle size (mm) vol: amount of water sampled per depth (l) TotalParticles: The number of particles seen in that volume (#) depth: The depth we are sampling (m)

We process the data slightly

particle_data_processed <- particle_data %>% 
  mutate(binsize = ub-lb, # size of particle bins
         nparticles = TotalParticles/vol, # particles normalized to volume (#/L)
         n_nparticles = nparticles/binsize, # particles normalized to volume and bin size (#/L/mm)
  )
particle_data_processed

I model the relationship between particles and depth as a power law function. That is the log of the particle size is linearly related to the log of the size and volume normalized particle numbers. To account for zeros, I use a poisson glm

myGlm <- function(df){
  glm(TotalParticles ~ log(lb), offset = log(vol * binsize), family = poisson, data = df)
}
particle_icp_psd <- particle_data_processed %>% nest(-depth, -profile) %>% 
  mutate(model = map(data, myGlm)) %>%
  mutate(tidied = map(model, tidy)) %>%
  unnest(tidied) %>%
  select(depth, profile, term, estimate) %>%
  spread(key = term, value = estimate) %>%
  rename(icp = `(Intercept)`, psd = `log(lb)`)

All elements of `...` must be named.
Did you want `data = c(lb, ub, vol, TotalParticles, binsize, nparticles, n_nparticles)`?

particle_icp_psd

One more thing. I also care about the particle flux. Flux is the sum of all of the particles times their sinking speed, times their mass. I relate flux to mass as follows Flux{in a bin of diamter D} = C_f * D ^ ag Cf and and ag are constants Cf = 4, ag = 0.23 And Total Flux = sum_D{Flux(D)} Flux is calculated based on the particles normalized to volume

#C_f = 4
#ag = 0.23
particle_flux <- particle_data_processed %>% 
  mutate(flux = (C_f * nparticles ^ ag)) %>%
  group_by(depth, profile) %>%
  summarize(Flux = sum(flux))
particle_flux

I combine the flux data with the PSD and intercept data

particle_ipf <- left_join(particle_icp_psd, particle_flux, by = "depth")
particle_ipf

pFlux <- ggplot(particle_ipf, aes(x = Flux, y = depth)) + scale_y_reverse() + geom_point()
pPSD <- ggplot(particle_ipf, aes(x = psd, y = depth)) + scale_y_reverse() + geom_point()
picp <- ggplot(particle_ipf, aes(x = icp, y = depth)) + scale_y_reverse() + geom_point()
plot_grid(pFlux, pPSD, picp, nrow = 1)

Now, for “reasons”, I want to smooth these profiles out and interpolate some new spectra. First, I use gams to model the three profiles.

gamPSD <- gam(psd ~ s(depth), data = particle_ipf)
gamicp <- gam(icp ~ s(depth), data = particle_ipf)
gamFlux <- gam(Flux ~ s(depth), family = "Gamma", data = particle_ipf)

predData <- tibble(depth = seq(from = 0, to = 1000, by = 25))
predPSD <- predict(gamPSD, predData, type = "response")
predicp <- predict(gamicp, predData, type = "response")
predFlux <- predict(gamFlux, predData, type = "response")

predDf <- tibble(depth = predData$depth, psd = predPSD, icp = predicp, Flux = predFlux)
predDf

pFlux <- ggplot(predDf, aes(x = Flux, y = depth)) + scale_y_reverse() + geom_point() + scale_x_continuous(limits = c(0, 200))
pPSD <- ggplot(predDf, aes(x = psd, y = depth)) + scale_y_reverse() + geom_point()
pInt <- ggplot(predDf, aes(x = icp, y = depth)) + scale_y_reverse() + geom_point() 
plot_grid(pFlux, pPSD, pInt, nrow = 1)

Removed 3 rows containing missing values (geom_point).

So next, I’m going to recreate spectra from psd and int. I will see if these recreated spectra approximate the observed “Flux” values.

Recall that psd and int relate to binsize and volume normalized particle numbers.

lb_vec = particle_data %>% pull(lb) %>% unique # the particle sizes
binsize_vec = particle_data_processed %>% pull(binsize) %>% unique # the particle sizes

make_spectrum <- function(icp, psd, bins = lb_vec, binsizes = binsize_vec){
  C_n = exp(icp)
  nnp = (C_n * bins ^ psd) * (binsizes)
  nnp
}

make_spectrum <- function(icp, psd, bins = lb_vec, binsizes = binsize_vec){
  nnp = exp(log(bins) * psd + icp)
  np = nnp * binsizes
  tibble(lb = bins, np)
}


predDf <- predDf %>% mutate(spec = map2(icp, psd, make_spectrum))

predDf[["spec"]][[5]] %>% mutate(flux = (C_f * np ^ag)) %>% summarize(Flux = sum(flux))

Now, I re-calculate flux from spec.

sumflux <- function(df){
  df %>% pull(flux) %>% sum
}
predDf2 <- predDf %>% 
  mutate(spec = map(spec, 
                    . %>% mutate(flux = (C_f * np ^ ag))
                    )) %>%
  #mutate(Flux = map(spec, ~ . %>% summarise(Flux = sum(flux))))
  mutate(Flux = map_dbl(spec, sumflux))

pPredFlux <- ggplot(data = predDf2, aes(y = depth, x = Flux)) + geom_point() + scale_y_reverse() + scale_x_continuous(limits = c(0, 200))
plot_grid(pFlux, pPredFlux)

Removed 3 rows containing missing values (geom_point).Removed 4 rows containing missing values (geom_point).

Unnest PredDf

predDf3 <- predDf2 %>% unnest(spec)

New names:
* `111` -> `111...11`
* `112` -> `112...12`
* `113` -> `113...13`
* `114` -> `114...14`
* `115` -> `115...15`
* ...

Disagg part

Data

synthetic_data <- predDf2 %>% 
  mutate(spec_only = map(spec, ~pull(., np)),
    prev_spec = lag(spec_only), 
         prev_Flux = lag(Flux),
         DF = prev_Flux - Flux,
         DFP = 1-DF/prev_Flux
    # I could have calulated DFP = Flux / prev_Flux, which is equivalent to this way.
  )

Variabiles

# mass of a 1mm particle
m1mm = 3.3 * 10^-6; #%g % Alldgedge 1998 % mass of 1mm particle
w1mm = 2; #% m/day # Alldredge and Gotschalk, methinks % sinking speed of 1mm particle
micron = 1e-6;
# fractle dimension
#C_f = 4 # Usual Way
#ag = 0.26

# So we can conform to usual expectations of alpha and gamma being positive
C_f = C_f # defined earlier
ag = ag # defined earlier

Cm = m1mm
Cw = w1mm
m_vec =  Cm * lb_vec ^ alpha;
w_vec = Cw * lb_vec ^ gamma;

test_abun_in <- synthetic_data %>% filter(depth == 25) %>% pull(spec) %>%.[[1]] %>% pull(np)

f_vec0 <- test_abun_in * m_vec * w_vec
F0 <- sum(f_vec0)

test_abun_out <- synthetic_data %>% filter(depth == 50) %>% pull(spec) %>% .[[1]] %>% pull(np)
f_vec1 <- test_abun_out * m_vec * w_vec
F1 <- sum(f_vec1)
DFP <- F1/F0
DFP

[1] 0.4871608

little_lb <- lb_vec[1] - (lb_vec[2] - lb_vec[1])/2 # size of the particle that the UVP can't see anymore. Eg, things actually shrink to this size but then they vanish from the UVP's view. Just leting it be like, the difference in size of the smallest two bins smaller than the smallest bin.

remin_shuffle <- function(abun_in, DFpct, DeltaZ = 10, size = lb_vec, Cm = m1mm, Cw = w1mm, lbv = lb_vec, mv = m_vec, wv = w_vec,
                          alpha = 0.52, gamma = 0.26, llb = little_lb){
 rn = abun_in * lb_vec
 ran = abun_in * lb_vec ^ alpha
 srn = sum(rn)
 sran = sum(ran)
 
 omega = lb_vec[1] ^ (2 * alpha)  * abun_in[1]/(lb_vec[1] ^ alpha - llb ^ alpha)
 #omega = lb_vec[1] ^ (2 * alpha)  * abun_in[1]/(llb ^ alpha - lb_vec[1] ^ alpha) # possible correction
 
 nmw = abun_in * mv * wv
 F1 = sum(nmw)
 DeltaF = (F1 * DFpct) - F1 # should be negative
 
 #Cr = DeltaF/ (Cm*(1+gamma/alpha) * DeltaZ * (sran + omega));
 Cr = DeltaF/ (Cm*  DeltaZ * ((1+gamma/alpha) * sran + omega)); # Possible correction
 #Cr = DeltaF/ (Cm*(1+gamma/alpha) * DeltaZ * (sran)); # Test
 CrCw = Cr/Cw
 
 phi = CrCw * (1+gamma/alpha) * ran * 
   DeltaZ/
   (c(llb, lbv)[1:length(lbv)]^alpha - lbv ^ alpha) # added extra parentheses
 
 Delta_nj_out = phi/lbv ^ gamma
 Delta_nj_in = c(phi[2:length(phi)],0)/c(lbv[2:length(phi)],1) ^ gamma
 #Delta_nj_in = -c(phi[2:length(phi)],0)/lbv ^ gamma
 Delta_nj_net = Delta_nj_in - Delta_nj_out # positive because out is negative
   
  return(list(Cr = Cr, phi = phi, dnet = Delta_nj_net, din =  Delta_nj_in, dout = Delta_nj_out))
}

Smooth Shuffle

remin_shuffle_spec <- function(abun_in, ...){
  core <- remin_shuffle(abun_in, ...)
  abun_in + core$dnet
}

dds_test <- 0.99
trs_spec <- remin_shuffle_spec(test_abun_in, dds_test, 2, llb = .08) # .076, .102, 0.08
trs_spec

 [1] 4.977118e+01 3.220295e+01 2.080654e+01 1.324086e+01 8.429664e+00 5.259372e+00 3.412537e+00 2.158349e+00 1.376862e+00 8.570189e-01
[11] 5.706750e-01 3.591209e-01 2.222774e-01 1.423506e-01 9.117745e-02 5.843491e-02 3.667519e-02 2.348329e-02 1.496326e-02 9.434999e-03
[21] 6.126706e-03 3.875805e-03 2.409464e-03 1.576325e-03 8.755244e-04

(pre <- sum(C_f * test_abun_in ^ ag))

[1] 763.788

(expected <- sum(C_f * test_abun_in ^ ag) * dds_test)

[1] 756.1501

(modeled <- sum(C_f * trs_spec ^ ag))

[1] 755.6298

(expected - pre)/pre

[1] -0.01

(modeled - pre)/pre

[1] -0.01068125

(modeled - expected)/expected

[1] -0.0006881335

llb seems to be pretty relevant for getting flux right – presumably I am not handling it correctly. I’d like things to adjust so that it is considered correctly

test_abun_in * dds_test

 [1] 4.995311e+01 3.238547e+01 2.091009e+01 1.329027e+01 8.457959e+00 5.266899e+00 3.419150e+00 2.159717e+00 1.377626e+00 8.557674e-01
[11] 5.705486e-01 3.586372e-01 2.217339e-01 1.419959e-01 9.090876e-02 5.825105e-02 3.653090e-02 2.338785e-02 1.489793e-02 9.387706e-03
[21] 6.096721e-03 3.855430e-03 2.395040e-03 1.568487e-03 8.751845e-04

sum((C_f * test_abun_in ^ ag)) * dds_test

[1] 756.1501

explore_llb <- tibble(llb = c(0, 0.01, 0.07, 0.08, 0.09, .1, .101, .1019, .10199)) %>%
  mutate(spec = map(llb, ~remin_shuffle_spec(test_abun_in, dds_test, 2, llb = .)))
  
explore_llb_unwrap <- explore_llb %>%
  mutate(spec = map(spec, ~data_frame(lb_vec = lb_vec, spec = .))) %>%
  unnest(spec) %>%
  pivot_wider(names_from = lb_vec, values_from = spec)

`data_frame()` is deprecated as of tibble 1.1.0.
Please use `tibble()` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

explore_llb_flux <- explore_llb %>%
  mutate(flux = map_dbl(spec, ~sum(C_f * . ^ ag))) %>%
  select(-spec)

explore_llb_2 <- left_join(explore_llb_flux, explore_llb_unwrap)

Joining, by = "llb"

  #
explore_llb_2

As we approach llb = lb_vec[1]

dds_test <- 0.99
trs_spec <- remin_shuffle_spec(test_abun_in, dds_test, 2, llb = .08) # .076, .102
trs_spec

 [1] 4.977118e+01 3.220295e+01 2.080654e+01 1.324086e+01 8.429664e+00 5.259372e+00 3.412537e+00 2.158349e+00 1.376862e+00 8.570189e-01
[11] 5.706750e-01 3.591209e-01 2.222774e-01 1.423506e-01 9.117745e-02 5.843491e-02 3.667519e-02 2.348329e-02 1.496326e-02 9.434999e-03
[21] 6.126706e-03 3.875805e-03 2.409464e-03 1.576325e-03 8.755244e-04

sum(C_f * test_abun_in ^ ag)

[1] 763.788

sum(C_f * trs_spec ^ ag)

[1] 755.6298

remin_smooth_shuffle <- function(abun_in, DFpct, Ipct = 0.9999, ...){
  # DFpct: Fractional mass retained between depths
  # Ipct: Fractional mass retained between iterations
  # ...: Passed to remin_shuffle
  IMirror <- 2  - Ipct
  
  abun_est = abun_in
  Fpct = DFpct # gets overwritten if we are iterating
  
  # If we are loossing flux, and we loose loose more flux than Ipct
  # Iterate remin shuffle only keeping ipct each time
  if(DFpct < Ipct){
    iters = floor(log(DFpct)/log(Ipct)) # why is this a ratio of log transformed values? 
    # I need to understand this before I can address the case where DFpct > IMirror
    iterFlux = Ipct^iters
    Fpct = 1-(iterFlux-DFpct)
    
    for (i in 1:iters){
      abun_est = remin_shuffle_spec(abun_in = abun_est, Ipct, ...)
    }
  }
  
  # If we are gaining flux, and we gain more than IMirror, iterate
  if(DFpct > IMirror){
    iters = floor(log(DFpct)/log(IMirror))
    iterFlux = IMirror^iters
    Fpct = 1 - (iterFlux - DFpct)
    
    for (i in 1:iters){
      abun_est = remin_shuffle_spec(abun_in = abun_est, IMirror, ...)
    }
  }
  
  # Deal with remainder. In the case where the loss is less than ipct, or greater than 2-ipct (Imirror), just do this part
  abun_est = remin_shuffle_spec(abun_in = abun_est, Fpct, ...)
  
  abun_est
}

DFP = 0.9
trss0 <- remin_shuffle_spec(DFpct = DFP , abun_in = test_abun_in) # .076, .102
trss0

 [1] 3.573079e+01 2.919141e+01 1.894655e+01 1.215562e+01 7.757633e+00 4.900530e+00 3.169375e+00 2.021353e+00 1.290119e+00 8.133358e-01
[11] 5.373678e-01 3.405726e-01 2.122546e-01 1.359706e-01 8.733901e-02 5.604447e-02 3.534732e-02 2.265133e-02 1.446004e-02 9.154131e-03
[21] 5.939993e-03 3.766087e-03 2.351742e-03 1.529026e-03 8.252953e-04

trss <- remin_smooth_shuffle(DFpct =  DFP, abun_in = test_abun_in) # .076, .102
trss

 [1] 3.713991e+01 2.901460e+01 1.883700e+01 1.208027e+01 7.715583e+00 4.869549e+00 3.151720e+00 2.009593e+00 1.283620e+00 8.086840e-01
[11] 5.345717e-01 3.389169e-01 2.112405e-01 1.353414e-01 8.694556e-02 5.580600e-02 3.520412e-02 2.256208e-02 1.440645e-02 9.121834e-03
[21] 5.919266e-03 3.753961e-03 2.344373e-03 1.522550e-03 8.206877e-04

# Confirm whether flux is treated appropriately
sum(C_f * test_abun_in ^ ag) * DFP

[1] 687.4092

sum(C_f * trss0 ^ ag)

[1] 670.6552

sum(C_f * trss ^ ag)

[1] 673.6107

tibble(test_abun_in, trss0, trss, lb_vec) %>% ggplot(aes(x = lb_vec)) + geom_point(aes(y = trss0)) + geom_point(aes(y = trss), shape = 2) + scale_y_log10() + scale_x_log10() + geom_point(aes(y = test_abun_in), shape = 1)

Graveyard for now

# I want to minimize this function here
shuffle_check <- function(abun_in, DFpct, C_f2 = C_f, ag2 = ag, ...){
  abun_out <- remin_smooth_shuffle(abun_in, DFpct, ...)
  flux_in <- sum(C_f2 * abun_in ^ ag2)
  flux_out <- sum(C_f2 * abun_out ^ ag2)
  DFPct_actual = flux_out/flux_in
  rmse <- (DFpct - DFPct_actual)^2
  rmse
}

find_DFP <- function(abun_in, ...){
  abun_in_loc = abun_in
  
  sc_wrap <- function(x, abun_in_loc){
    shuffle_check(abun_in = abun_in_loc, DFpct = x)
  }
  
  adjDfp <- optimize(sc_wrap, c(0, 2), abun_in_loc = test_abun_in)
  adjDfp
}


# remin_smooth_shuffle_fix <- function(abun_in, DFpct, ...){
#   
# }

DFP = 0.8
shuffle_check(DFpct = DFP , abun_in = test_abun_in)

[1] 0.001137407

find_DFP(abun_in = test_abun_in)

$minimum
[1] 0.9999975

$objective
[1] 2.269221e-13

Data

First, lets calculate Dfpct for each pair of depths, and lets also stagger spec

dummy_function <- function(DFpct, spec, ...){
  DFpct * spec
}
dummy_function(0.8, test_abun_in)

 [1] 4.036615e+01 2.617008e+01 1.689704e+01 1.073961e+01 6.834714e+00 4.256080e+00 2.762950e+00 1.745226e+00 1.113233e+00 6.915292e-01
[11] 4.610493e-01 2.898079e-01 1.791789e-01 1.147442e-01 7.346163e-02 4.707156e-02 2.951992e-02 1.889927e-02 1.203873e-02 7.586025e-03
[21] 4.926643e-03 3.115499e-03 1.935386e-03 1.267464e-03 7.072198e-04

synthetic_remineralized <- synthetic_data %>% .[-1,] %>% mutate(pred_spec = map2(prev_spec, DFP, remin_smooth_shuffle))
synthetic_remineralized

Did it work?

tibble(
actual = synthetic_remineralized[["spec_only"]][[4]],
predicted = synthetic_remineralized[["pred_spec"]][[4]],
previous = synthetic_remineralized[["prev_spec"]][[4]],
lb = lb_vec
) %>% gather(key = "var", value = "np", actual, predicted, previous) %>%
  ggplot(aes(x = lb, y = np, shape = var, col = var), ) + geom_point(size = 3, alpha = 0.75) + scale_x_log10() + scale_y_log10() + labs(x = "# Particles", y = "Size (mm)", col = "Model", shape = "Model")

Frigging finally. Consistant everything I think.

Looks promising! Lets be more systematic about this next! synthetic_remineralized_con Ok. If I calculate flux assuming alpha = 2.3 and gamma = 1.3, I get this stronger than traditional arch effect. I guess because flux attenuates faster than it should.

If I say that alpha + gamma = 0.63 as best fit flux, then I loose more big particles than small ones. More resetting instnaces and now things look really strange.

As of april 21 2020, I’m using alldredge alpha and gamma

But also my different flux profiles are no longer adding up.

As above, but this time, with adjusted DFP to account for “leakage”

Leaks

I’m having a leakage problem. I get nearly, but not quite the right amount of flux out. What if I address this by running remin smooth shuffle and then I optimize the DFPct that it gets so that it returns the right flux output.

I ought to actually figure out what is wrong, but this should give a close answer for most analyses.

Optimization function: What I want it to do. Takes a structure and a DFP. Varies the DFP sent to remin smooth shuffle so the acutal DFP loss is what we want

We need a function where we can pass the DFPct that we are using and the actual DFPct

shuffle_tune <- function(DFpct_toRemin, abun_in,  DFpct_target,...){
  abun_out <- remin_smooth_shuffle(abun_in, DFpct_toRemin, ...)
  flux_in <- sum(C_f_global * abun_in ^ ag_global)
  flux_out <- sum(C_f_global * abun_out ^ ag_global)
  DFPct_actually_happened <- flux_out/flux_in
  rmse <- (DFPct_actually_happened - DFpct_target)^2
  rmse
}

stOpt <- optimize(shuffle_tune, c(0, 2), abun_in = test_abun_in, DFpct_target = 0.9)
stOpt

$minimum
[1] 0.9154147

$objective
[1] 3.926557e-11

So for a DFpct of 0.9, I should actually feed it 0.915

Apply stopt over the data set.

optFun <- function(abun_in, DFpct){
  opt <- optimize(shuffle_tune, c(0, 2), abun_in = abun_in, DFpct_target = DFpct)
  opt$minimum
}

Test optimization

plot(
  seq(from = 0.5, to = 1, by = 0.05),
  map_dbl(seq(from = 0.5, to = 1, by = 0.05), ~optFun(testSpectrum, .))
  )

Apply optimization

synthetic_data_fixDfp <- synthetic_data %>% filter(depth > 0) %>% mutate(use_this_DFpct = map2_dbl(prev_spec, DFP, optFun))

synthetic_data_fixDfp %>% ggplot(aes(x = DFP, y = use_this_DFpct)) + geom_point() + geom_abline(v = 1,h = 1)

Ignoring unknown parameters: v, h

synthetic_data_fixDfp %>% ggplot(aes(x = DFP, y = use_this_DFpct - DFP)) + geom_point() + geom_abline(v = 1,h = 1)

Ignoring unknown parameters: v, h

Redo remineralization

synthetic_remineralized <- synthetic_data_fixDfp %>% .[-1,] %>% mutate(pred_spec = map2(prev_spec, use_this_DFpct, remin_smooth_shuffle))
synthetic_remineralized

Analyzing remineralized data

Indicating which regions are actually loosing flux the whole time. Lets first look at particle numbers, observed minus expected

dfpPlt <- synthetic_remineralized %>% ggplot(aes(y = depth, x = DFP, col = DFP < 1)) + geom_point() + scale_y_reverse() 
fluxPlt <- synthetic_remineralized %>% ggplot(aes(y = depth, x = Flux, col = DFP < 1)) + geom_point() + scale_y_reverse()
plot_grid(dfpPlt, fluxPlt)

synthetic_remineralized_proc <- synthetic_remineralized %>% mutate(TP = map_dbl(spec, sum), pred_TP = map_dbl(pred_spec, sum), dif_TP = TP-pred_TP) 

dtpP <- synthetic_remineralized_proc %>% ggplot(aes(y = depth, x = dif_TP, col = DFP < 1)) + geom_point() + scale_y_reverse()
ptpP <- synthetic_remineralized_proc %>% ggplot(aes(y = depth, x = pred_TP, col = DFP < 1)) + geom_point() + scale_y_reverse() + scale_x_log10()
tpP <- synthetic_remineralized_proc %>% ggplot(aes(y = depth, x = TP, col = DFP < 1)) + geom_point() + scale_y_reverse() + scale_x_log10()
plot_grid(tpP, ptpP, dtpP)

Now that I’ve fixed this for positive values, we predict absurd numbers of particles much of the time. Lets just look at decreasing parts.

synthetic_remineralized_proc_decreasing <- synthetic_remineralized_proc %>% filter(DFP <= 1)

dtpP <- synthetic_remineralized_proc_decreasing %>% ggplot(aes(y = depth, x = dif_TP, col = DFP < 1)) + geom_point() + scale_y_reverse()
ptpP <- synthetic_remineralized_proc_decreasing %>% ggplot(aes(y = depth, x = pred_TP, col = DFP < 1)) + geom_point() + scale_y_reverse() + scale_x_log10()
tpP <- synthetic_remineralized_proc_decreasing %>% ggplot(aes(y = depth, x = TP, col = DFP < 1)) + geom_point() + scale_y_reverse() + scale_x_log10()
fractpP <- synthetic_remineralized_proc_decreasing %>% ggplot(aes(y = depth, x = dif_TP/TP, col = DFP < 1)) + geom_point() + scale_y_reverse()
plot_grid(tpP, ptpP, dtpP, fractpP)

Oh dear. The model predicts profound particle number attenuation. Actually there isn’t so everything actually just follows flux.

Maybe there’s some better metric. Like fraction of flux on particles smaller than 53 micron or something.

This may have been a dead end. I need to reflect again.

useRow <- which(synthetic_remineralized$DFP > 1)[1]

tibble(
actual = synthetic_remineralized[["spec_only"]][[useRow]],
predicted = synthetic_remineralized[["pred_spec"]][[useRow]],
previous = synthetic_remineralized[["prev_spec"]][[useRow]],
lb = lb_vec
) %>% gather(key = "var", value = "np", actual, predicted, previous) %>%
  ggplot(aes(x = lb, y = np, shape = var, col = var), ) + geom_point(size = 3, alpha = 0.75) + scale_x_log10() + scale_y_log10()

Ah. Somethings off here.

test_abun_in

 [1] 5.045769e+01 3.271260e+01 2.112130e+01 1.342451e+01 8.543393e+00 5.320100e+00 3.453687e+00 2.181533e+00 1.391541e+00 8.644115e-01
[11] 5.763117e-01 3.622598e-01 2.239736e-01 1.434302e-01 9.182704e-02 5.883945e-02 3.689990e-02 2.362409e-02 1.504842e-02 9.482531e-03
[21] 6.158304e-03 3.894373e-03 2.419233e-03 1.584330e-03 8.840247e-04

trss00 <- remin_smooth_shuffle(DFpct =  1, abun_in = test_abun_in) # .076, .102
trss00

 [1] 5.045769e+01 3.271260e+01 2.112130e+01 1.342451e+01 8.543393e+00 5.320100e+00 3.453687e+00 2.181533e+00 1.391541e+00 8.644115e-01
[11] 5.763117e-01 3.622598e-01 2.239736e-01 1.434302e-01 9.182704e-02 5.883945e-02 3.689990e-02 2.362409e-02 1.504842e-02 9.482531e-03
[21] 6.158304e-03 3.894373e-03 2.419233e-03 1.584330e-03 8.840247e-04

trss00 <- remin_shuffle_spec(DFpct =  1.01, abun_in = test_abun_in) # .076, .102
trss00

 [1] 5.193038e+01 3.306471e+01 2.133878e+01 1.355140e+01 8.621969e+00 5.362057e+00 3.482118e+00 2.197551e+00 1.401683e+00 8.695191e-01
[11] 5.802061e-01 3.644286e-01 2.251455e-01 1.441762e-01 9.227584e-02 5.911895e-02 3.705516e-02 2.372136e-02 1.510726e-02 9.515371e-03
[21] 6.180135e-03 3.907202e-03 2.425982e-03 1.589860e-03 8.898977e-04

13 April 2020

Hmm. Because of the power law nature of particles and size essentially “all” of the particles are created at each depth. Essentially particle number isn’t a great proxie for disaggregation.

I have a couple of ideas.

One is to look at <53 miron vs >53 micron particles. One could ask: How much biomass or flux moves from one side of that line to the other?

I’d also like to get at things like how particle production relates to flux attenuation. Like how much of flux attenuation is different with the small particles than if there weren’t these small particles.

I could imagine calculating C_r everywhere and then just running a purely prism model propigating downward with that Cr. Seeing what flux is at some key threshold, and then comparing that to actual flux. I’d expect that essentially all of the attenuation woudl be through the production and removal of small particles.

In the back of my mind though, my assumptions about small particles are bugging me. What if they are more recalcitrant than big particles. Then all of the models sort of don’t work, right?

20 April 2020

Some unnesting

synthetic_remineralized_concise <- synthetic_remineralized %>%
  mutate(spec2 = map2(spec, prev_spec, ~tibble(.x, prev_np = .y))) %>%
  mutate(spec2 = map2(spec2, pred_spec, ~tibble(.x, pred_np = .y))) %>%
  select(-c(spec, prev_spec, pred_spec, spec_only))

synthetic_remineralized_unnested <- synthetic_remineralized_concise %>% unnest(spec2)

I wonder what is up with these name reassignments?

Goal: Estimate attenuation of large and small particles. Compare that to total attenuation. (Hopefully they sum to the same thing). Do this for both the model, and for the observed data.

Calculate flux transfer to small particles from disaggregation which should be

ObservedSmallFlux - PredictedSmallFlux Which should equal: PredictedBigFlux - ObservedFlux

sr02 <- synthetic_remineralized_unnested %>%
  mutate(prev_flux = C_f * prev_np ^ ag,
         pred_flux = C_f * pred_np ^ ag)

sr02Tot <- sr02 %>%
  group_by(depth) %>%
  summarize(DF = first(DF), DFP = first(DFP),
            Flux = sum(flux), Prev_Flux = sum(prev_flux), Pred_Flux = sum(pred_flux))

sr02Small <- sr02 %>% 
  filter(lb <= 0.53) %>%
  group_by(depth) %>%
  summarize(Flux = sum(flux), Prev_Flux = sum(prev_flux), Pred_Flux = sum(pred_flux))

sr02Big <- sr02 %>% 
  filter(lb > 0.53) %>%
  group_by(depth) %>%
  summarize(Flux = sum(flux), Prev_Flux = sum(prev_flux), Pred_Flux = sum(pred_flux))

sr02All <- sr02Tot %>%
  left_join(sr02Small, by = "depth", suffix = c("", "_Small")) %>%
  left_join(sr02Big, by = "depth", suffix = c("", "_Big")) %>%
  mutate(osps = Flux_Small - Pred_Flux_Small, obpb =  Pred_Flux_Big - Flux_Big) # And of course they are not equal
head(sr02All)

Ok. Osps ~= obpb, which means that big flux lost about equals small flux gained due to disaggregation (or whatever proces makes things deviate from the model)

sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = osps, col = DFP)) + geom_point() + scale_y_reverse()

Huh. Look at that. Sometimes osps goes negative, which is to say there is apparent aggregation there. I think.

sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = DFP, x = osps)) + geom_point() + scale_y_reverse()

Ok. What do I actually want to know here. I think osps devided by the flux attenuation of the big particles (predicted or observed). I want to know what fraction of big particles become small particles vs just dissolve.

Or maybe something else too?

sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = DF, x = osps)) + geom_point() + scale_y_reverse()

ospsPlt <- sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = osps/DF)) + geom_point() + scale_y_reverse() + labs(x = "Disaggregation") + theme_cowplot() + geom_vline(xintercept = 0, col = "gray50")
ospsPlt

Huh.Look at this. We have disaggregation happening in the surface, at ~500m (after an uptick in flux, which gets erased, so maybe a zooplankton thing) and again at the bottom of the OMZ. We have net aggregation in the OMZ core, for the most part.

This is cool!

plot_grid(fluxPlt + labs(col = "IsFluxDecreasing?") + theme_cowplot() + theme(legend.position = "none"),
          ospsPlt + theme_cowplot())

I think the stuff below is misguided and the thing above, especially combined wiht illustrations of where flux increases, will be sufficiently informative.

Disaggregation are positive values here, btw. Aggregation are negative values.

sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = osps/(Prev_Flux_Big- Flux_Big))) + geom_point() + scale_y_reverse() + scale_x_log10()

Why is the above ever more than one? More stuff disaggregates than disapears? If it disapears, doesn’t it disaggregate?

plot_prev_big <- sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = Prev_Flux_Big-Flux_Big)) + geom_point() + scale_y_reverse() 
plot_big <- sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = Prev_Flux_Big - Pred_Flux_Big)) + geom_point() + scale_y_reverse() 
plot_pred_big <- sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = Pred_Flux_Big)) + geom_point() + scale_y_reverse() + scale_x_log10()

cowplot::plot_grid(plot_prev_big, plot_big, plot_pred_big)

Huh. Sometimes the big flux goes up, even when total flux goes down, both in observed and predicted space. I guess this makes sense if we have attenuation, but a flattening of the curve.

sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = osps/(Prev_Flux_Big- Pred_Flux_Big))) + geom_point() + scale_y_reverse() + scale_x_log10()

All this and more tomorrow.

plot_grid(dfpPlt, fluxPlt)

p <- sr02All %>% filter(DFP <= 1) %>% ggplot(aes(y = depth, x = osps/(Prev_Flux_Big - Flux_Big))) + geom_point() + scale_y_reverse()
plotly::ggplotly(p)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

875 m is the most clea example. Lets look at the profiles and see what is happening there.

tibble( actual = synthetic_remineralized[[“spec_only”]][[4]], predicted = synthetic_remineralized[[“pred_spec”]][[4]], previous = synthetic_remineralized[[“prev_spec”]][[4]], lb = lb_vec ) %>% gather(key = “var”, value = “np”, actual, predicted, previous) %>% ggplot(aes(x = lb, y = np, shape = var, col = var), ) + geom_point(size = 3, alpha = 0.75) + scale_x_log10() + scale_y_log10()

synthetic_remineralized_unnested %>% filter(depth == 875, lb < 0.2) %>% select(lb, actual = np, previous = prev_np, predicted = pred_np) %>%
gather(key = "var", value = "np", actual, previous, predicted) %>%
   ggplot(aes(x = lb, y = np, shape = var, col = var), ) + geom_point(size = 3, alpha = 0.75) + scale_x_log10() + scale_y_log10()

synthetic_remineralized_unnested %>% filter(depth == 875, lb > 10) %>% select(lb, actual = np, previous = prev_np, predicted = pred_np) %>%
gather(key = "var", value = "np", actual, previous, predicted) %>%
   ggplot(aes(x = lb, y = np, shape = var, col = var), ) + geom_point(size = 3, alpha = 0.75) + scale_x_log10() + scale_y_log10()

synthetic_remineralized_unnested %>% filter(depth == 875) %>% mutate(previous = prev_np - np, predicted = pred_np - np) %>% select(lb, previous, predicted) %>%
gather(key = "var", value = "np",previous, predicted) %>%
   ggplot(aes(x = lb, y = 1/(np), shape = var, col = var), ) + geom_point(size = 3, alpha = 0.75) + scale_x_log10()

Huh. The curve appears to be flattening more than we would expect from random chance. But there are still more particles because of shenanigans in the smallist bin.

So these positive numbers might in fact be places with negative osps (but also flattening and so an opposite of usual change in large particles)

I wonder how even to address this.